
R version 4.1.2 (2021-11-01) -- "Bird Hippie"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> ## Code to replicate Table 1 of Abrajano, Elmendorf, and Quinn. "Measuring
> ## Perceived Skin Color: Spillover Effects and Likert-type Scales". JOP.
> 
> 
> library(clubSandwich)
Registered S3 method overwritten by 'clubSandwich':
  method    from    
  bread.mlm sandwich
> 
> 
> sink(file="AbrajanoElmendorfQuinn-Table1-output.txt")
> 
> ## Part 1
> ## E[y | d = W] > E[y | d = B]
> ##
> ## where the expectation is over all respondents assigned to M1 or M2,
> ## all photos, and all target positions
> 
> mydata <- read.csv("ScaleRaceSpring2017clean.csv")
> 
> ## just black distractor photos
> mydata.b <- mydata[mydata$condition %in% c("P1B", "P2B"),]
> ## just white distractor photos
> mydata.w <- mydata[mydata$condition %in% c("P1W", "P2W"),]
> 
> ## keep must the target photos
> target.inds <- c(6, 7, 10, 12, 15, 17, 19, 20, 22, 23)
> keep.vars <- c(paste("X", target.inds, "_Q339", sep=""),
+                paste("X", target.inds, "_Q329", sep=""))
> 
> mydata.b <- mydata.b[, keep.vars]
> mydata.w <- mydata.w[, keep.vars]
> 
> ## reshape into long format
> mydata.b.long <- reshape(mydata.b, direction="long", varying=1:20, v.names="MM")
> mydata.w.long <- reshape(mydata.w, direction="long", varying=1:20, v.names="MM")
> 
> ## keep only obs with non-missing values for MM
> mydata.b.long <- na.omit(mydata.b.long)
> mydata.w.long <- na.omit(mydata.w.long)
> 
> lm.b.out <- lm(MM ~ 1, data=mydata.b.long)
> lm.w.out <- lm(MM ~ 1, data=mydata.w.long)
> 
> tab.b <- coef_test(lm.b.out, vcov="CR2", cluster=mydata.b.long$id,
+                    test="Satterthwaite")
> 
> tab.w <- coef_test(lm.w.out, vcov="CR2", cluster=mydata.w.long$id,
+                    test="Satterthwaite")
> WminusB.diff <- tab.w[1,1] - tab.b[1,1]
> WminusB.var <- tab.w[1,2]^2 + tab.b[1,2]^2
> WminusB.SE <- sqrt(WminusB.var)
> WminusB.z <- WminusB.diff / WminusB.SE
> WminusB.pval <- 2*(1-pnorm(abs(WminusB.z)))
> 
> cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
> cat("Part 1:  E[y | d = W] - E[y | d = B] > 0\n")
> cat("\n")
> cat("----------------------------------------------------------------\n")
> cat("Avg. MM score among respondents with white distractors\n")
> print(as.data.frame(tab.w), digits=4)
> cat("\nAvg. MM score among respondents with black distractors\n")
> print(as.data.frame(tab.b), digits=4)
> cat("\n\n----------------------------------------------------------------\n")
> cat("WminusB.diff     WminusB.SE     WminusB.z    p-value (2-tailed)\n")      
> cat("----------------------------------------------------------------\n")
> cat(WminusB.diff, "       ", WminusB.SE, "    ", WminusB.z, "       ",
+     WminusB.pval, "\n")
> cat("----------------------------------------------------------------\n")
> cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> ## Part 2
> ## E[y_s | d = W] > E[y_s | d = B]
> ## where the expectation is over all respondents assigned to M1 or M2,
> ## all photos, and conditional on s = 6.
> 
> mydata <- read.csv("ScaleRaceSpring2017clean.csv")
> 
> ## just black distractor photos
> mydata.b <- mydata[mydata$condition %in% c("P1B", "P2B"),]
> ## just white distractor photos
> mydata.w <- mydata[mydata$condition %in% c("P1W", "P2W"),]
> 
> ## keep just the target photos
> target.inds <- c(6)
> keep.vars <- c(paste("X", target.inds, "_Q339", sep=""),
+                paste("X", target.inds, "_Q329", sep=""))
> 
> mydata.b <- mydata.b[, keep.vars]
> mydata.w <- mydata.w[, keep.vars]
> 
> ## reshape into long format
> mydata.b.long <- reshape(mydata.b, direction="long", varying=1:2, v.names="MM")
> mydata.w.long <- reshape(mydata.w, direction="long", varying=1:2, v.names="MM")
> 
> ## keep only obs with non-missing values for MM
> mydata.b.long <- na.omit(mydata.b.long)
> mydata.w.long <- na.omit(mydata.w.long)
> 
> lm.b.out <- lm(MM ~ 1, data=mydata.b.long)
> lm.w.out <- lm(MM ~ 1, data=mydata.w.long)
> 
> tab.b <- coef_test(lm.b.out, vcov="CR2", cluster=mydata.b.long$id,
+                    test="Satterthwaite")
> 
> tab.w <- coef_test(lm.w.out, vcov="CR2", cluster=mydata.w.long$id,
+                    test="Satterthwaite")
> WminusB.diff <- tab.w[1,1] - tab.b[1,1]
> WminusB.var <- tab.w[1,2]^2 + tab.b[1,2]^2
> WminusB.SE <- sqrt(WminusB.var)
> WminusB.z <- WminusB.diff / WminusB.SE
> WminusB.pval <- 2*(1-pnorm(abs(WminusB.z)))
> 
> cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
> cat("Part 2:  E[y_6 | d = W] - E[y_6 | d = B] > 0\n")
> cat("\n")
> cat("----------------------------------------------------------------\n")
> cat("Avg. MM score among respondents with white distractors\n")
> print(as.data.frame(tab.w), digits=4)
> cat("\nAvg. MM score among respondents with black distractors\n")
> print(as.data.frame(tab.b), digits=4)
> cat("\n\n----------------------------------------------------------------\n")
> cat("WminusB.diff     WminusB.SE     WminusB.z    p-value (2-tailed)\n")      
> cat("----------------------------------------------------------------\n")
> cat(WminusB.diff, "       ", WminusB.SE, "    ", WminusB.z, "       ",
+     WminusB.pval, "\n")
> cat("----------------------------------------------------------------\n")
> cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")
> 
> 
> 
> 
> 
> 
> 
> 
> ## Part 3
> ## E [y | d = W, i in I_M1 ]  -  E [y | d = B, i in I_M1 ] >
> ##  E[y | d = W, i in I_M2 ]  -  E [y | d = B, i in I_M2 ] 
> ## where the expectations are over all respondents assigned to M1 or M2,
> ## all photos, and all target positions.
> 
> mydata <- read.csv("ScaleRaceSpring2017clean.csv")
> 
> ## just black distractor photos
> mydata.b <- mydata[mydata$condition %in% c("P1B", "P2B"),]
> ## just white distractor photos
> mydata.w <- mydata[mydata$condition %in% c("P1W", "P2W"),]
> 
> ## keep just the target photos
> target.inds <- c(6, 7, 10, 12, 15, 17, 19, 20, 22, 23)
> keep.vars.template <- paste("X", target.inds, "_Q329", sep="")
> keep.vars.notemplate <- paste("X", target.inds, "_Q339", sep="")
> 
> mydata.b.m1 <- mydata.b[, keep.vars.template]
> mydata.b.m2 <- mydata.b[, keep.vars.notemplate]
> mydata.w.m1 <- mydata.w[, keep.vars.template]
> mydata.w.m2 <- mydata.w[, keep.vars.notemplate]
> 
> ## reshape into long format
> mydata.b.m1.long <- reshape(mydata.b.m1, direction="long",
+                             varying=1:10, v.names="MM")
> mydata.b.m2.long <- reshape(mydata.b.m2, direction="long",
+                             varying=1:10, v.names="MM")
> mydata.w.m1.long <- reshape(mydata.w.m1, direction="long",
+                             varying=1:10, v.names="MM")
> mydata.w.m2.long <- reshape(mydata.w.m2, direction="long",
+                             varying=1:10, v.names="MM")
> 
> ## keep only obs with non-missing values for MM
> mydata.b.m1.long <- na.omit(mydata.b.m1.long)
> mydata.b.m2.long <- na.omit(mydata.b.m2.long)
> mydata.w.m1.long <- na.omit(mydata.w.m1.long)
> mydata.w.m2.long <- na.omit(mydata.w.m2.long)
> 
> lm.b.m1.out <- lm(MM ~ 1, data=mydata.b.m1.long)
> lm.b.m2.out <- lm(MM ~ 1, data=mydata.b.m2.long)
> lm.w.m1.out <- lm(MM ~ 1, data=mydata.w.m1.long)
> lm.w.m2.out <- lm(MM ~ 1, data=mydata.w.m2.long)
> 
> tab.b.m1 <- coef_test(lm.b.m1.out, vcov="CR2", cluster=mydata.b.m1.long$id,
+                    test="Satterthwaite")
> tab.b.m2 <- coef_test(lm.b.m2.out, vcov="CR2", cluster=mydata.b.m2.long$id,
+                    test="Satterthwaite")
> 
> tab.w.m1 <- coef_test(lm.w.m1.out, vcov="CR2", cluster=mydata.w.m1.long$id,
+                    test="Satterthwaite")
> tab.w.m2 <- coef_test(lm.w.m2.out, vcov="CR2", cluster=mydata.w.m2.long$id,
+                    test="Satterthwaite")
> 
> WminusB.m1.diff <- tab.w.m1[1,1] - tab.b.m1[1,1]
> WminusB.m1.var <- tab.w.m1[1,2]^2 + tab.b.m1[1,2]^2
> WminusB.m1.SE <- sqrt(WminusB.m1.var)
> WminusB.m1.z <- WminusB.m1.diff / WminusB.m1.SE
> WminusB.m1.pval <- 2*(1-pnorm(abs(WminusB.m1.z)))
> 
> WminusB.m2.diff <- tab.w.m2[1,1] - tab.b.m2[1,1]
> WminusB.m2.var <- tab.w.m2[1,2]^2 + tab.b.m2[1,2]^2
> WminusB.m2.SE <- sqrt(WminusB.m2.var)
> WminusB.m2.z <- WminusB.m2.diff / WminusB.m2.SE
> WminusB.m2.pval <- 2*(1-pnorm(abs(WminusB.m2.z)))
> 
> WminusB.diff.diff <- WminusB.m1.diff - WminusB.m2.diff
> WminusB.diff.diff.var <- WminusB.m1.var + WminusB.m2.var
> WminusB.diff.diff.SE <- sqrt(WminusB.diff.diff.var)
> WminusB.diff.diff.z <- WminusB.diff.diff / WminusB.diff.diff.SE
> WminusB.diff.diff.pval <- 2*(1-pnorm(abs(WminusB.diff.diff.z)))
> 
> 
> cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
> cat("Part 3: (E [y | d = W, i in I_M1 ]  -  E [y | d = B, i in I_M1 ]) -\n")
> cat("         (E[y | d = W, i in I_M2 ]  -  E [y | d = B, i in I_M2 ]) > 0\n") 
> cat("\n")
> cat("----------------------------------------------------------------\n")
> cat("Avg. MM score among respondents with white distractors & no template\n")
> print(as.data.frame(tab.w.m1), digits=4)
> cat("\nAvg. MM score among respondents with black distractors & no template\n")
> print(as.data.frame(tab.b.m1), digits=4)
> cat("\n\n")
> cat("Avg. MM score among respondents with white distractors & template\n")
> print(as.data.frame(tab.w.m2), digits=4)
> cat("\nAvg. MM score among respondents with black distractors & template\n")
> print(as.data.frame(tab.b.m2), digits=4)
> 
> cat("\n\n----------------------------------------------------------------\n")
> cat("WminusB.m1.diff  WminusB.m1.SE  WminusB.m1.z   p-value (2-tailed)\n")      
> cat("----------------------------------------------------------------\n")
> cat(WminusB.m1.diff, "       ", WminusB.m1.SE, "    ", WminusB.m1.z, "       ",
+     WminusB.m1.pval, "\n")
> cat("----------------------------------------------------------------\n")
> 
> cat("\n\n----------------------------------------------------------------\n")
> cat("WminusB.m2.diff  WminusB.m2.SE  WminusB.m2.z   p-value (2-tailed)\n")      
> cat("----------------------------------------------------------------\n")
> cat(WminusB.m2.diff, "       ", WminusB.m2.SE, "    ", WminusB.m2.z, "       ",
+     WminusB.m2.pval, "\n")
> cat("----------------------------------------------------------------\n")
> 
> cat("\n\n****************************************************************\n")
> 
> cat("----------------------------------------------------------------\n")
> cat("WminusB.dif.dif  WminusB.dif.dif  WminusB.dif.dif.z   p-value (2-tailed)\n")      
> cat("----------------------------------------------------------------\n")
> cat(WminusB.diff.diff, "       ", WminusB.diff.diff.SE, "    ", WminusB.diff.diff.z, "       ",
+     WminusB.diff.diff.pval, "\n")
> cat("----------------------------------------------------------------\n")
> 
> cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")
> 
> 
> 
> sink()
> 
> proc.time()
   user  system elapsed 
  1.442   0.127   1.734 
